An exploration of nighttime lights and the Dark-Sky Movement in the US
Electric lights have benefited society immensely since their invention in 1879. Modern work, recreation, and transportation are made possible largely due to artificial lighting. However, when artificial outdoor lighting becomes inefficient and unnecessary, it is known as light pollution, with negative impacts on both wildlife and humans [1].
Concerns about the impacts of light pollution on observational astronomy led to the 1988 formation of the nonprofit now known as DarkSky International [2]. Since then, DarkSky has grown into a recognized global authority on light pollution with a mission to safeguard communities and wildlife against light pollution [3]. The organization’s International Dark Sky Places program allows communities, parks, and protected areas to seek certification for efforts to preserve and protect their dark skies [4]. While many communities implement outdoor lighting policies without seeking recognition, a DarkSky certification is an indicator of places that have made dark sky conservation a priority.
This project aims to explore trends in nighttime light levels first from a national contiguous US perspective, and through more detailed case studies of individual cities. Three DarkSky certified communities [4] with different sizes and designation dates have been selected:
These three communities will each be paired with a community of similar size that has fewer, if any, ordinances controlling light pollution. Each community pair will be examined to see if there are observable differences in nighttime light emissions.
The following data were used.
library(tidyverse)
library(leaflet)
library(kableExtra)
library(htmlwidgets)
library(widgetframe)
knitr::opts_chunk$set(widgetframe_widgets_dir = 'widgets' )
knitr::opts_chunk$set(cache=TRUE) # cache the results for quick compiling
library(blackmarbler)
library(sf)
library(terra)
library(tigris)
library(tidyterra)
library(tidycensus)
library(patchwork)
library(piggyback)
library(ggiraph)A US Census TIGER shapefile was downloaded to provide geographic state boundaries. A separate script was used to download and save annual nighttime lights data for the contiguous US (‘black_marble_data_download.R’). Raster data was downloaded for plotting purposes, and extracted mean nighttime lights values were also downloaded for analysis. Nighttime lights are radiance values with units of \(\frac{n W}{cm^2 sr}\).
# Define national spatial extent
nonconus <- c("Guam", "Hawaii", "Alaska",
"Commonwealth of the Northern Mariana Islands",
"United States Virgin Islands", "American Samoa", "Puerto Rico")
states_sf <- states() |>
filter(!NAME %in% nonconus) |>
select(NAME, geometry) |>
st_transform(crs = "epsg:4326")
# Dissolve individual state boundaries to get CONUS outline
conus <- st_union(states_sf)|>
st_as_sf()
## Read in national raster data for 2024 (to be used for leaflet map later)
pb_download("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t2024.tif",
repo = "cstaebell/final-project-cstaebell")
us_2024_rast <- rast("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t2024.tif") |>
mask(conus)
# create log transformed version for plotting
log_us_2024 <- app(us_2024_rast, fun = log1p)
#---
# Download and process annual US mean NTL data
#---
years = 2014:2024
conus_means = numeric(length(years))
for (i in 1:length(years)) {
pb_download(paste("CONUS_VNP46A4_NearNadir_Composite_Snow_Free_qflag_mean_t",
years[i], ".Rds", sep = ""), repo = "cstaebell/final-project-cstaebell")
annual_mean <- readRDS(paste("CONUS_VNP46A4_NearNadir_Composite_Snow_Free_qflag_mean_t",
years[i], ".Rds", sep = "")) |>
pull(ntl_mean)
conus_means[i] <- annual_mean
}
conus_df <- data.frame(year = years, mean_ntl = conus_means)
#---
# Download and process mean NTL data for states
#---
state_means <- matrix(nrow = 49, ncol = length(years))
# Read in annual data for states and create a matrix of state means
for (i in 1:length(years)) {
pb_download(paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_mean_t",
years[i], ".Rds", sep = ""), repo = "cstaebell/final-project-cstaebell")
annual_mean <- readRDS(paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_mean_t",
years[i], ".Rds", sep = ""))
state_means[,i] <- annual_mean[,2]
# check that state names are ordered the same way for each year
if (i == 1) {
state_names <- annual_mean[,1]
} else {
if (sum(state_names != annual_mean[,1]) != 0) {
print("Different state names")
}
}
}
# Create data frame of state means and add percent change column
state_means_df <- data.frame(state_means)
colnames(state_means_df) <- paste("ntl_", years, sep = "")
state_means_df <- state_means_df |>
mutate(state = state_names,
pct_change = (ntl_2024 - ntl_2014) / ntl_2014 * 100)
# Combine state NTL data with geographic data
states_ntl <- states_sf |>
inner_join(state_means_df, by = c("NAME" = "state"))For case study analyses, population data and local geographies for each Dark Sky Place were downloaded from the US Census.
# Download population data
# US population to search for companion cities
us_pop <- get_acs(geography = "place", state = NULL,
variables = "B01003_001", key = census_api_key)
# State populations for the states in which the selected DarksKy communities reside
az_pop <- get_acs(geography = "place", state = "AZ",
variables = "B01003_001", key = census_api_key)
il_pop <- get_acs(geography = "place", state = "IL",
variables = "B01003_001", key = census_api_key)
tx_pop <- get_acs(geography = "place", state = "TX",
variables = "B01003_001", key = census_api_key)
# Download local geographies: selected Dark Sky Places
flagstaff <- places(state = "AZ", cb = TRUE, year = 2024) |>
filter(NAME == "Flagstaff") |>
st_transform(crs = "epsg:4326") |>
left_join(us_pop, by = "GEOID")
homerglen <- places(state = "IL", cb = TRUE, year = 2024) |>
filter(NAME == "Homer Glen") |>
st_transform(crs = "epsg:4326") |>
left_join(us_pop, by = "GEOID")
beecave <- places(state = "TX", cb = TRUE, year = 2024) |>
filter(NAME == "Bee Cave") |>
st_transform(crs = "epsg:4326") |>
left_join(us_pop, by = "GEOID")Candidate cities to pair with each selected Dark Sky community were found by filtering US ACS population data to return a listing of municipalities with populations within 1% of each Dark Sky community population. Final comparison cities were selected based on population match and state. Preference was given to cities within states that are not known to have laws regarding light pollution or dark sky protection, as compiled by Schultz [12].
# Get listing of suitable comparison cities
flagstaff_sim_pop <- us_pop |> filter(estimate >= (flagstaff$estimate - 0.005*flagstaff$estimate)
& estimate <= (flagstaff$estimate + 0.005*flagstaff$estimate)) |>
separate(NAME, into = c("city", "state"), sep = ", ")
homerglen_sim_pop <- us_pop |> filter(estimate >= (homerglen$estimate - 0.005*homerglen$estimate)
& estimate <= (homerglen$estimate + 0.005*homerglen$estimate)) |>
separate(NAME, into = c("city", "state"), sep = ", ")
beecave_sim_pop <- us_pop |> filter(estimate >= (beecave$estimate - 0.005*beecave$estimate)
& estimate <= (beecave$estimate + 0.005*beecave$estimate)) |>
separate(NAME, into = c("city", "state"), sep = ", ")
# Download/process comparison city data
# Flagstaff counterpart
scranton <- places(state = "PA", cb = TRUE, year = 2024) |>
filter(NAME == "Scranton") |>
left_join(us_pop, by = "GEOID") |>
st_transform(crs = "epsg:4326")
# Homer Glen counterpart
belvidere <- places(state = "IL", cb = TRUE, year = 2024) |>
filter(NAME == "Belvidere") |>
left_join(us_pop, by = "GEOID") |>
st_transform(crs = "epsg:4326")
# Bee Cave counterpart
clanton <- places(state = "AL", cb = TRUE, year = 2024) |>
filter(NAME == "Clanton") |>
left_join(us_pop, by = "GEOID") |>
st_transform(crs = "epsg:4326")The final city pairs are as follows:
| Pair | Community | Population |
|---|---|---|
| 1 | Flagstaff, AZ | 76,333 |
| 1 | Scranton, PA | 76,074 |
| 2 | Homer Glen, IL | 24,516 |
| 2 | Belvidere, IL | 24,510 |
| 3 | Bee Cave, TX | 8,861 |
| 3 | Clanton, AL | 8,862 |
The annual US nighttime lights rasters were then subsetted to create a raster and dataframe of mean values for each city.
# Subset rasters to each area
cities <- list(flagstaff, homerglen, beecave, scranton, belvidere, clanton)
city_rasters <- vector("list", length = length(cities))
# Load files for each year
for (i in 1:length(years)){
pb_download(paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t",
years[i], ".tif", sep = ""), repo = "cstaebell/final-project-cstaebell")
rast_file <- paste("US_VNP46A4_NearNadir_Composite_Snow_Free_qflag_t",
years[i], ".tif", sep = "")
year_rast <- rast(rast_file)
# Loop over cities
for (j in 1:length(cities)) {
# Crop and mask to city extent
current_rast <- year_rast |>
terra::crop(cities[[j]]) |>
mask(cities[[j]])
# Create list entry during first iteration, then append in subsequent iterations
if (i == 1) {
city_rasters[[j]] <- list(current_rast)
} else {
city_rasters[[j]][[length(city_rasters[[j]]) + 1]] <- current_rast
}
}
}
# Stack city lists
city_rasters <- lapply(city_rasters, rast)
# Create separate rasters for each city
city_var_names <- paste0(c("flagstaff", "homerglen", "beecave", "scranton",
"belvidere", "clanton"), "_rast")
for (c in 1:length(city_var_names)) {
assign(city_var_names[c], city_rasters[[c]])
}
# Calculate annual means
annual_means <- vector("list", length = length(cities))
for (i in 1:length(city_rasters)) {
ntl_vals <- data.frame(values(city_rasters[[i]])) |>
summarize(across(t2014:t2024, ~ mean(.x, na.rm = TRUE)))
annual_means[[i]] <- ntl_vals
}
# Create mean dataframes for each city
city_df_var_names <- paste0(c("flagstaff", "homerglen", "beecave", "scranton",
"belvidere", "clanton"), "_means")
for (i in 1:length(annual_means)) {
df <- annual_means[[i]] |>
pivot_longer(cols = t2014:t2024) |>
mutate(name = years) |>
rename(year = name, mean_ntl = value)
assign(city_df_var_names[i], df)
}However, simply looking at annual radiance means for each city may be more of a proxy for overall development or population density than for lighting policy. Land use/land cover data can be used to help make comparisons across communities more meaningful.
Although land cover changes over the period from 2014 to 2024, for the purposes of this project, only one year of land cover data was considered for each community. Data availability is one reason for this simplification – annual data is not available for all study areas for all years, and the MODIS science team urges caution when using land cover layers for 2021 and beyond due to lack of updates to the classification training database [11]. Therefore, 2019 data was used for all cities except for Clanton, which used 2020 data.
Land cover rasters were downloaded and cropped to each city. The urban areas in each city were identified, and then annual mean nighttime lights values were calculated for that urban area.
# Download land cover data from GitHub Release
pb_download("MCD12Q1.061_LC_Type1_doy2019001000000_aid0001.tif",
repo = "cstaebell/final-project-cstaebell")
pb_download("MCD12Q1.061_LC_Type1_doy2020001000000_aid0001.tif",
repo = "cstaebell/final-project-cstaebell")
land_use_rast_2019 <- rast("MCD12Q1.061_LC_Type1_doy2019001000000_aid0001.tif")
land_use_rast_2020 <- rast("MCD12Q1.061_LC_Type1_doy2020001000000_aid0001.tif")
# Subset land cover rasters to each city
city_ntl_rasters <- list(flagstaff_rast, homerglen_rast, beecave_rast, scranton_rast,
belvidere_rast, clanton_rast)
land_use_rasters <- vector("list", length = length(cities))
urban_ntl_dfs <- vector("list", length = length(cities))
for (i in 1:length(cities)) {
# Specify which data to use (2019 for all cities but Clanton)
if (i == 6) {
land_use_rast <- land_use_rast_2020
} else {
land_use_rast <- land_use_rast_2019
}
# Crop and mask to city extent
current_rast <- land_use_rast |>
terra::crop(cities[[i]]) |>
mask(cities[[i]])
# Save to list for future use
land_use_rasters[[i]] <- current_rast
# Mask all but urban land use (land use code = 13)
urban_area <- current_rast
urban_area[urban_area != 13] <- NA
city_urban <- city_ntl_rasters[[i]] |>
mask(urban_area)
# Calculate mean NTL for urban areas
urban_ntl_dfs[[i]] <- data.frame(values(city_urban)) |>
summarize(across(t2014:t2024, ~ mean(.x, na.rm = TRUE))) |>
pivot_longer(cols = t2014:t2024) |>
mutate(name = years) |>
rename(year = name, mean_urban_ntl = value)
}
# Create separate dataframes from list object
urban_ntl_var_names <- paste0(c("flagstaff", "homerglen", "beecave", "scranton",
"belvidere", "clanton"), "_urban_ntl")
for (i in 1:length(urban_ntl_dfs)) {
assign(urban_ntl_var_names[i], urban_ntl_dfs[[i]])
}The 2024 nighttime lights across the contiguous US are visualized below.
ggplot() +
geom_spatraster(data = log_us_2024) +
scale_fill_gradient2(low = "black",
mid = "yellow",
high = "red",
midpoint = 4.5,
na.value = "transparent") +
labs(title = "Nighttime Lights: 2024") +
coord_sf() +
theme_void() +
theme(legend.position = "none")
Mean nighttime light radiance for the entire contiguous US is plotted below, where a generally increasing trend can be observed over the past ten years.
ggplot(conus_df) +
geom_point(aes(x = year, y = mean_ntl)) +
geom_smooth(aes(x = year, y = mean_ntl)) +
labs(x = "Year", y = expression("Mean Radiance (nW cm"^-2~"sr"^-1*")"),
title = "Annual US Mean Radiance (2014 - 2024)") +
theme_light()
State means for 2024 are plotted in the figure below. Annual radiance tends to be higher for eastern states, which reflects known patterns of population density.
# Create ggplot map for ggiraph
st_2024_mean_map <- ggplot() +
geom_sf_interactive(data = filter(states_ntl, NAME != "District of Columbia"),
aes(fill = ntl_2024,
tooltip = paste(NAME, "\n", round(ntl_2024, 2)),
data_id = NAME)) +
theme_void() +
theme(legend.position = "bottom",
legend.direction = "horizontal") +
scale_fill_viridis_c() +
labs(title = "Mean Annual Nighttime Lights by State, 2024",
fill = expression("Mean Radiance (nW cm"^-2~"sr"^-1*")"),
caption = "Data: NASA Black Marble & US Census Bureau")
# Transform into interactive ggiraph object
int_st2024mean <- girafe(ggobj = st_2024_mean_map,
options = list(opts_hover(css = "stroke:cyan;stroke-width:2px"),
opts_hover_inv(css = "opacity:0.8;"),
opts_tooltip(
css = "background-color:white;color:black;padding:15px;
border-radius:10px;
font-family:Arial;font-size:15px;"
)))
int_st2024mean